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ABSTRACT 

Using iV-body simulations with a large set of massless test particles we compare the 
predictions of two theories of violent relaxation, the well known Lynden-Bell theory 
and the more recent theory by Nakamura. We derive "weakened" versions of both 
theories in which we use the whole equilibrium coarse-grained distribution function fi 
as a constraint instead of the total energy constraint. We use these weakened theo- 
ries to construct expressions for the conditional probability i"Q( r ) that a test particle 
initially at the phase-space coordinate r would end-up in the i'th macro-cell at equi- 
librium. We show that the logarithm of the ratio Rij{r) = K^t)/ Kj(t) is directly 
proportional to the initial phase-space density /o( r ) f or the Lynden-Bell theory and 
inversely proportional to fo (r) for the Nakamura theory. We then measure Rij (t) us- 
ing a set of iV-body simulations of a system undergoing a gravitational collapse to 
check the validity of the two theories of violent relaxation. We find that both theories 
are at odds with the numerical results, both qualitatively and quantitatively. 

Key words: methods: numerical - galaxies: haloes - galaxies: statistics - galaxies: 
kinematics and dynamics 



1 INTRODUCTION 

A statistical theory that successfully describes the process 
of a collisionless gravitational collapse has been a longstand- 
ing open problem. The ultimate goal of such a theory is to 
predict the (semi-) equilibrium state of a gravitational sys- 
tem that is described by the Vlasov equation, given its ini- 
tial phase-space density /o(r, v). The implications of such 
a theory to cosmology and astrophysics are immense. For 
example, such a theory may be able to explain the origin of 
cusps found in simulations of CDM haloes, the universality 
of their shape, and possibly connect them to the initial cold 
phase-space density distribution. 

The first attempt to construct such a theory was 
made by Lynden-B ell in his pioneering work from 1967 
llLvnden-Bell lll967l hereafter the LB67 theory) which ac- 
tually coined the phrase "violent relaxation" to describe the 
relaxation of such systems. However, already in that paper 
Lynden-Bell predicted that his theory will not be applicable 
to large parts of the system due to the "incompleteness of 
violent relaxation" . The fluctuations of the gravitational po- 
tential, which drive the system towards an equilibrium state 
fade out too soon for the system to explore the full config- 
uration space. Therefore the final equilibrium state picked 
by the system is not necessarily the one that maximises the 
entropy. The discrepancy becomes larger as we go from the 
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centre of the system to its outer parts where the gravita- 
tional fluctuations are much smaller. Accordingly, Lynden- 
Bell assumed that the complete relaxation is limited to a 
sphere of radius Ri around the system's centre of mass. 

Soon after this paper, many A-body simulations 
were run to check its validity. Most of these simula- 
tions were of one-dimensional models of either plane 



sheets or of spherical shells JCohen & Lecar 1 


196c; 


Henonlll968: CuDcrman at el. 


RM 


Goldstein et al. 1 


1969; 


Lecar & Cohen Ill97ll: iTanekusa Ill987f). In most cases the 



LB67 theory was only partially correct at best. A common 
outcome of the "water-bucket" initial conditions (in which 
there is only one level of phase-space density surrounded 
by a vacuum) was the "core-halo" structure in which the 
halo took most of the energy of the system, leaving the core 
(which contained most of the mass) degenerate. Compared 
with predictions, the phase-space density was too high at 
very low energies, too low at intermediate energies, and os- 
cillated at high energies. One can find a good fit for the core 
using the degenerate solution of the LB67 theory, and a rea- 
sonable fit for the halo using the non-degenerate solution - 
but then it is difficult to find a clean way to decide which 
particles should be described by which one of the two solu- 
tions. Consequently, the theory loses much of its predictive 
power. This ambiguity was attributed to the incompleteness 
of the relaxation: the inner parts where the relaxation was 
effective could be described by a LB67 solution, unlike the 
outer parts which were only partially relaxed. 

The disagreement between the predictions of the LB67 
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theory and experiments, in addition to some disturbing con- 
ceptual issues it possesses (such as its infinite mass pre- 
diction in 3D and phase-space densities segregation - see 
Sec. 12.31 for further detail s ), hav e led people to cons i der a l- 
ternatives to it. e.g..lShu I ll978h:IStiavelli fc Bertin I Jl987t): 
ISpergel fc Hernauist Hl992h:lKuH et al. I l)l997l) : lNakamural 
feOOCft : iTrenti fc Bertin I J20051) . In all these theories the 
equilibrium state is assumed to be the most probable state 
of the system, which is found by maximising the entropy 
under the appropriate constraints. The difference between 
these theories is mainly due to the definition of entropy they 
use and the constraints under which it is maximised. 

For the purpose of this paper we roughly divide these 
theories into two groups: in the first group we have theories 
which are based on a more fundamental approach to the 
problem, thereby incorporating (almost) all of the dynami- 
cal constraints when maximising the entropy. In this group 
we ess entially have two t heories: the t heory of lLy ndcn-Bell 
1 1967) and the theories o f lKull et al. lll997l) aiid lN^Jtainunu 
i 2000), which as we shall see, are basically the same theory. 



In the second group we have theories that use more heuris- 
tic , ad-hoc like, cons t raints and entropy definitions, suc h 
as IStiavelli fc Bertin I jl987l): ISpergel fc Hernauist I Jl992>) : 
ITrenti fc Bertin I (120051) . or lHansen et al. I (120051) wh o use 
non-ex tensive statistical mechanics, and to some extent lshu I 
(1978). Non of these theories, for example, take into ac- 
count the conservation of the phase-space volume by the 
dynamics, as required by Liouville's theorem. As a rule of 
thumb, the predictions of the second group are in a better 
agreement with simulations and observations than the first 
group. However, this is achieved at the expense of increasing 
the number of assumptions and free parameters in the the- 
ory. In this paper we are interested in examining how well 
the fundamental assumptions of statistical mechanics apply 
in violent relaxation, we shall therefore concentrate on the 
first group of the ories: the LB67 theory and the theory of 
iNakamura I <l2000l) (hereafter NK00). 

The purpose of this paper is three-fold: firstly, we wish 
to compare the predictions of the LB67 and NK00 theories 
with numerical experiments done in 3D. To our knowledge 
this has not yet been done directly, in particular when it 
comes to the NK00 theory. While the ID case shares many 
common characteristics with the 3D case, it is known that 
dimensionality may play an important role in systems un- 
dergoing violent relaxation. For example, the 3D system has 
more degrees of freedom, and may therefore experience a 
more efficient mixing. Indirect support for the NK 00 the- 
ory in 3D is found in iMerrall fc Henrikseiil ll2003l) . Here 
the authors show that the velocity distribution function of 
isolated systems after violent relaxation is well fitted by a 
single Gaussian, as is predicted by the NK00 theory for the 
non-degenerate case. Nevertheless, we would like to perform 
a more direct test of the predictions of this theory. 

Secondly, acknowledging the fact that these theories are 
largely incorrect, we would still like to see if certain aspects 
of them are true. In other words, we would like to see if the 
approach of maximising the entropy under the proper con- 
straints is at all useful in the weakest possible sense for a dis- 
sipationless self-gravitating system. Thirdly, as we shall see 
in Sec. |21 the NK00 theory and the LB67 theory have a very 
similar structure, incorporating the same set of constraints, 
while using fundamentally different entropy definitions. We 



would therefore like to know which definition describes vio- 
lent relaxation better. We feel that somehow the question of 
how entropy should be defined is more fundamental than the 
full statistical theory around it. For example, a local defini- 
tion of entropy can be used in a dynami cal theory t hat d e- 
scribes the approach to equilibrium [e.g. IChavanis I 1^998)]. 

The structure of this paper is as follows: in Sec.[5]we re- 
view the LB67 and NK00 theories and derive them using the 
information-theory approach. We introduce the conditional 
and joint probabilities that describe the path of a test par- 
ticle and relate them to the different definitions of entropy 
between the two theories. We then list some of the prob- 
lems and questions with these theories and re-derive "weak- 
ened" versions of the theories. We show how the predictions 
of these theories can be tested in a numerical experiment 
that to a large extent avoids many of the pitfalls we would 
have encountered had we tried to test the LB67 and NK00 
theories in their original form. In Sec. we describe how 
the conditional probabilities can be measured in an iV-body 
simulation using a large number of test particles. Then in 
Sec. 0] we describe the numerical simulations that we have 
performed to measure these probabilities and in Sec. we 
present the results of the simulations. Our conclusions are 
presented in Sec.|S] 



2 THEORETICAL BACKGROUND - TWO 
THEORIES FOR VIOLENT RELAXATION 

In this section we will focus on the LB67 theory and a more 
recent alternative to it - the NK00 theory. The aim of these 
theories is to predict the final equilibrium state of a colli- 
sionless gravitating system, given its initial state. The state 
of the system is completely determined by its phase-space 
density function (DF) /(r, v, £), which is governed by the 
Vlasov equation 



9tf + V • drf - V$ ■ d v f = 



(1) 



Here <l?(r, t) is the gravitational potential, calculated self- 
consistently from p(r, t) = J/(r, v, t)d 3 v via Poisson's equa- 
tion. 



2.1 Derivation of the NKOO theory 

The NKOO theory was originally derived in the framework 
of information theory. Here we shall repeat its derivation 
for comparison with the LB67 theory. We first formulate 
the violent relaxation process as an experiment with a well 
defined set of possible results Ui and a set of corresponding 
probabilities pi. The most probable result is then the one 
that maximises the Shannon entro py S = y~l, Pi log P i under 
some prescribed constraints on pi l|,Tavnes Ill957al lbf). 

Let /o(t) be the initial phase-space density of the sys- 
tem, with r = (r, v) being a phase-space coordinate. We 
toss a test particle into the initial phase-space according to 
the probability distribution po (r) which is defined by 



Po(t) 



M 



Mr) , 



(2) 



and let it move under gravity just like any phase-space ele- 
ment. We define p{r,t) to be the probability distribution of 
finding the test particle at time t > 0. The conservation of 
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phase-space volume guarantees that p(r, t) = f(r,t)/M for 
all t > 0. 

Next, we divide phase-space into macro-cells i — 
1, 2, 3, . . . of volume Co, and define the coarse-grained prob- 
ability pi as the probability of finding the test particle in 
the i'th macro-cell when the system reaches an equilibrium. 
From the above discussion it is clear that pi is equal to fi /M 
with fi being the coarse-grained DF in the macro-cell i at 
equilibrium. 

A possible result of our experiment is the pair (r, i) 
which specifies the initial location of the test particle and the 
number of the cell where it is finally found. With each such 
pair we associate the joint probability- distribution Pi(r). 
Pi{r)d a r is the probability that the particle was initially 
in a small patch of volume d 6 r around r, and ended up in 
the i'th macro-cell. The Shannon entropy of our experiment 
is then 



(3) 



Snk = — d e 'r Pi{r) log Pi{r) 



Maximising Snk gives us the most probable Pi(r), from 
which we can get 



fi = M Pi (T)d 6 T 



(4) 



Before maximising Snk, however, we must write down the 
constraints on the probabilities Pi(r). The first constraint 
stems from the initial conditions and from the fact that Pi(r) 
is a joint probability distribution: 



^Pi{-r) =Po(t) 



M 



Mr) 



(5) 



The second constraint comes from the conservation of phase- 
space volume under the collisionless dynamics: the overall 
initial phase-space volume of all phase-space elements that 
ended up in cell i will be exactly Q: 



Mr) 



(6) 



This constraint guarantees that the NK00 theory does not 
violate Liouville's theorem. Finally, from the conservation 
of energy, the total energy in fi must be equal to the initial 
energy of the system E 



E 



2 2 ^ Ir, - r 



(7) 



with (r», Vj) being the phase-space coordinates of the centre 
of the z'th cell. 

The initial conditions and the phase-space volume 
constraints © can be neatly written in terms of the condi- 
tional probability distribution Ki(r), defined by 



Ki{r) 



Pi(r) 



M 



(8) 



Poij) " /o(r) 

Ki(r) is the probability that a test particle which is initially 
at t will end-up at the i'th cell. Then the initial conditions 
and volume preservation constraints can be written as 



i 

Jd a rK x (r) 



= 1 



(9) 



(10) 



while the coarse-grained DF is given by 

fi = Jf {r)K l {r)d\ . (11) 

Using the well known technique of Lagrange multipliers it 
can now be easily shown that the Ki(r) for which 5 Snk = 
is 



Ki{ T ) = e 



3ej-«(-7-)-Aj//[j(-r) 



(12) 



Here a = vf/2 + $(ri) is the specific energy in the i'th cell 
and P,8{r),\i are the Lagrange multipliers of the energy 
constraint, the initial conditions constraint and the volume 
preservation constraint respectively. 



2.2 The LB67 theory: equal-mass discretisation 
vs. equal-volume discretisation 

The NK00 theory was originally derived within the infor- 
mation theory approach which was outlined in the previous 
section, while the LB67 theory was derived using a combi- 
natorial approach. Their different predictions, however, are 
not related to the different frameworks in which they were 
derived, but rathe r to the diffe r ent de finition of entropy they 
use. As noted bv iNakamura I ll200Cf) . the LB67 theory can 
be recovered in the framework of information theory if in- 
stead of defining the entropy in terms of the joint-probability 
Pi(r), one defines it in terms of the conditional probability 
Ki(r): 



S LB = -Y, (fr Ki(r) log Ki(r) 



(13) 



Indeed, maximising Slb with respect to the constraints JHJ- 
yields the expression 



Ki(r) = e-«»W'«-'W->' , 



(14) 



with f3, S(t), Xi being the Lagrange multipliers as in Eq. 112H . 
Using the volume preservation constraint (|SJ it is easy to 
express Xi in terms of the other unknowns, and plugging 
this into Eq. 1111 . after a trivial re-definition of 8 we obtain 



fi = 



/o>0 



-/3/o(T')e»-<5(V')d 6 T' 



(15) 



which is identical to the expression for the coarse-grained 
DF in the LB67 theory. 

On the other hand, lArad fc Lynden-BelTI J2005) have 
shown that the NK00 theory can be derived in a combi- 
natorial approach, where the entropy is defined simply (up 
to an additive constant) as the logarithm of the number of 
micro-states which comply with a given macro-state. While 
in the LB67 theory one discretises phase-space by consider- 
ing elements of equal volume, the NK00 theory is recovered 
if one considers elements of equal mass. Therefore we see 
that using the joint probability distribution in the informa- 
tion theory approach is equivalent to equal-mass discretisa- 
tion, whereas the use of conditional probability distribution 
is equivalent to equal-volume discretisation. The discretisa- 
tio n of phase-space in to cells of equal mass was already done 
bv lKull et al~T(ll997l). The NK00 theor y is therefore essen- 
tially the same as the lKull et al. I <ll997f) theory and for that 
reason we only consider one of them in this paper. 
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The above analogy can also be demonstrated by modi- 
fying the test particle experiment that was used in the pre- 
vious section to derive the NKOO theory. Instead of letting 
Po(t) - the probability distribution for finding the test par- 
ticle initially at r - be proportional to fo(r), we can devise 
an alternative experiment in which the particle has an equal 
probability of being everywhere in phase-space, provided, of 
course, that the overall phase-space volume V is finite. In 
other words, we let Po (t) = 1/V. As Pq (t) is a con- 
stant, the conditional probability k\ lb \t) is proportional 
to the joint-probability and therefore the usual Shannon en- 
tropy Eq. will yield the LB67 results. The connection to 
the equal-volume discretisation versus equal-mass discreti- 
sation difference that is found in the combinatorial approach 
is now evident: in the LB67 experiment the probability of 
initially finding the particle in a given patch of phase-space 
is proportional to the volume of the patch, whereas in the 
NKOO experiment it is proportional to the mass within that 
patch. 

Before concluding this section we note an important dif- 
ference between the resultant conditional probability in the 
LB67 theory and the NKOO theory. Whereas in the LB67 
theory the coupling between the initial coordinates r and 
the final coordinate i is in the /3/o(r)ej term in the expo- 
nent, in the NKOO theory it is in the Aj//o(r) term. This 
different coupling will be extensively discussed in the fol- 
lowing sections. 

2.3 Empirical and conceptual problems in the 
theories 

The statistical theory of violent relaxation (both LB67 and 
NKOO) suffers from several problems and open questions. In 
the next two sections we briefly describe some of these prob- 
lems and sketch a numerical experiment which will enable 
us to answer some of these questions. 

1. Incomplete relaxation Firstly, and perhaps most impor- 
tantly, the relax ation is never comp lete - indeed, as was al- 
ready noticed bv lLvnden-Bell I (Il967l) . after a few oscillations 
on the dynamical timescale of (Gp)~ 1 ^ 2 it is over. This does 
not give the system enough time to probe the configuration 
space, and as a result the system will settle down in a state 
that does not necessarily maximise the entropy. 
Nevertheless, it has often been argued that in the central 
regions of the collapsed object, where p is high and the dy- 
namical time scales are short, violent relaxation is efficient 
and we may therefore expect the system to approach the 
equilibrium solution in these regions. However, it is not clear 
how this claim can be verified quantitatively. The equilib- 
rium state depends on the gravitational potential, which in 
turn depends on the phase-space density also in the outer 
parts of the system, where the relaxation is incomplete. 

2. Definition of entropy As demonstrated in the last two 
sections, there exist (at least) two equally plausible ways of 
defining the entropy. These two definitions yield two differ- 
ent predictions for the equilibrium state, and therefore at 
least one of them is wrong. It would be interesting to check 
numerically which entropy (if either) better describes violent 
relaxation. 

3. Additional hidden constraints It is possible that in the 
violent relaxation process there exist a set of preserved quan- 



tities that are not taken into account in the theories when 
maximising the entropy. In such a case, the actual config- 
uration space of the system is much more limited and its 
maximal entropy state may be different from the calculated 
one. Moreover, this state would have a smaller entropy than 
the maximal entropy which is calculated without the con- 
straints. Such a mechanism may prove to be another reason 
for the incomplete violent relaxation. 

One such uncounted constraint is the total angular momen- 
tum of the system, which has not been taken into account 
in the present derivations of the LB67 and NKOO theo- 
ries. In principle, however, it can be taken into account by 
adding the appropriate Lagrange multipliers, with the cost 
of adding some additional complication to the final result. 
Another such invariant is 

(16) 

which was suggested bv lStiavelli fc Bertin I l)l987h to be ap- 
prox imately conserved upon a phenomenological basis [see 
also iTrenti fc Bertin I J2005)] . Further examples of invari- 
ants and their influen c e on the dynamic s can be found in 
iMoutarde et al I <1995l) : lHenriksen I (I2004T) . 
Finally the Hamiltonian nature of the phase-space flow gives 
rise to a set of constraints that are more difficult to han- 
dle. Such flows which can always be described by a canoni- 
cal transformation (r, v) — > (r', v'), necessarily preserve the 
so called Poincare integral invariants which are an infinite 
set of invariants. The circulation-like inte grals which were 
already mentioned in iLvnden-Bell I il967f) are a subset of 
these invariants. The inclusion of such invariants in a sta- 
tistical theory of violent relaxation is a mathematical chal- 
lenge, which to the best of our knowledge is yet to be over- 
come. 



2.4 The weakened versions of LB67 and NKOO 

To study these three points we have devised a numerical 
experiment in which the conditional probability Ki (r) is di- 
rectly measured. It is then compared to weakened versions 
of the LB67 and NKOO theories in which the whole equilib- 
rium coarse-grained DF f, is taken as a constraint, thereby 
replacing the total energy constraint. To see how this con- 
struction may shed light on the above points, let us first 
derive the weakened versions of the LB67 and NKOO theo- 
ries. 

We start with the NKOO theory, replacing the energy 
constraint £J with the "final conditions" Eq. 1111 . which is 
now taken as a constraint. The Lagrange multiplier j3 is thus 
replaced with a set of multipliers £j, and the functional we 
wish to maximise is 



Ink = -Y, <?T f {T)Ki{T)\n[f (T)Ki{T)] 
+ f\s(T)Ki(T) + XiKi(r) 
+ Zifo(T)Ki(T)\ dr . (17) 
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The functional derivative of I with respect to Ki(r) is then 
= -/ ( r ){ln[/ (r)^(r)] + l}+<5(r) 

+ Ai + /o(T)6 , 

and equating it to zero we obtain 

Intf,(r) = i£i+£ i + . Ai 



(18) 



Mr) SI fo(r) 



(19) 



after trivial redefinitions of S(t) and £j. 

To derive the equivalent result for the LB67 theory, we 
replace fo(r)Ki{r) — * Ki(r) in the definition of the entropy 
in Eq. p7f. while leaving the constraints intact. We thus 
obtain a new functional I lb- Differentiating it with respect 
to Ki(r) we get 



SIl 



lnKi(T) + l + 5(T) + \i + f (T)Zi , (20) 



«Q(t) 

and therefore the extremum solution is 

]nKi(T) = 6(T) + \i + f (T)£i 



(21) 



after redefining Ai — > Ai + 1. As noticed at the end of section 
12.21 the most striking difference between the two theories is 
the coupling between the initial coordinate r and the final 
coordinate i. This can be summarised as follows 



truth in either of them. If neither of the theories passes the 
test then there is an extremely small chance that the idea 
of maximising the entropy, at least as defined in the NK00 
and LB67 theories, has any physical relevance. On the other 
hand, if one of these theories passes the test it would allow 
us to decide what is the preferable form of the entropy. 

The proposed test also partly circumvents the difficul- 
ties in points 1 and 3: firstly, by picking i, j and r to be well 
inside the collapsed object, we are assured that the region of 
phase-space that we are measuring has undergone the maxi- 
mal mixing which the simulation provides. Of course it does 
not mean that this is the maximal mixing that is necessary 
for the entropy to reach its predicted maximum - but it is 
"as good as it can get". 

Secondly, as the above predictions for Rij (r) were made 
on the basis of theories which use the final fi as a constraint, 
any failure of the theories in passing the test could not be at- 
tributed to any unknown constraint which involve fi (such 
as energy and total angular momentum). In other words, 
the inclusion of such constraints would not change the pre- 
dictions of the theories for Rij(r). Other constraints, which 
involve Ki(r) itself, such as possibly the constraints due to 
the conservation of circulation integrals, may still affect the 
predictions. 



K[ NK \r) = A (NK \r)Bl NK) 
Kl LB) (r) = A (LB \r)B\ LB) e 



e /oM 

(LB 



Mr) 



(22) 
(23) 



We see very prominent differences between the two theo- 
ries, which may be detectable numerically if one measures 
Ki(r) directly. However, to find the theoretical predictions 
for A(r),Bi and d one has to solve the constraints Eqs. JO] 
I1UI 1111 . which is a non-trivial task. This can be largely 
avoided if we focus on the ratio 



R t] {r) = K t {T)/K 3 {r) 



(24) 



for some fixed i,j coordinates. The predictions of the NK00 
and LB67 theories for this quantity are: 



r\? k \t) = 



B 



(NK) 



R 



(is), 



r) = 



(NK) 

i 

(LB) 



(NK) (NK) 
c i j 

7oM 



(LB)_(LBY 



B 



(LB)' 



(25) 
(26) 



We notice that in both cases the only dependence of Rij(r) 
on r is through /o(t), and therefore by plotting Rij as a 
function of fo(r) we can see if either of the theories hold. 
Specifically, the LB67 theory predicts that In Rij (r) is linear 
in /o(t) whereas the NK00 theory predicts that it is linear 
in 1//o(t). This is the central idea of this paper. 

This is of course only a sufficient test. By passing it, 
we are not guaranteed that Ki(r) is given by Eq. I22H or 
Eq. 12311 . Moreover, even if this is the case, it is still not 
the form of the original NK00 and LB67 theories which use 
the much weaker energy constraint instead of the full coarse- 
grained fi constraint. But this is also the advantage of using 
such a test: we already know that these theories are wrong 
to a large extent, in particular at the outer parts of the 
collapsed object. Therefore just because this test is a very 
weak test, it allows us to see if there is a little grain of 



3 MEASURING CONDITIONAL 
PROBABILITIES IN AN TV-BODY 
SIMULATION 

To test the predictions of the LB67 and NK00 theories we 
ran a set of iV-body simulations of a gravitational collapse. 
The goal of these simulations was to measure the condi- 
tional probability Ki(r) from which the ratio Rij(r) = 
Ki(r) / Kj(r) can be calculated. In what follows we explain 
how this probability is measured. 



3.1 Using test particles to measure Ki(r) 

Theoretically, the conditional probabilities Ki(r) can be 
measured in an JV-body simulation using a large set of Nt 
test particles which are placed initially in a small patch of 
phase-space around r. The gravitational mass of the test 
particles is set to zero so that they do not affect the evolu- 
tion of the system - but only trace it. Then when the system 
relaxes one can estimate the phase-space density of the test 
particles by re-assigning to them a mass of m t = 1/Nt, such 
that their total mass is unity. The resultant density in the 
macro-cell i is then simply Ki(r). 

In practise, however, this procedure might fail as it is 
well known that in almost any iV-body simulation the tra- 
jectory of each particle exponentially departs from the ex- 
act traje ctory of the Vlasov equation due to chaos (see for 
example Irfeggie Ill99ll : iQuinlan fc Tremaine Ill992ft . These 
deviations, however, may only have a small effect on statis- 
tical measurements like the average density of some region 
in space, since the errors of individual particles tend to can- 
cel out each other. On the other hand, this cancellation, is 
not guaranteed when considering a group of test particles 
which are initially located in a tiny patch of phase-space. 
Such a configuration may be very sensitive to the particular 
realisation of the massive particles since all the trajectories 
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of the test particles are very close to each other and they can 
all be influenced simultaneously by a single massive particle. 
As a result, two simulations that use exactly the same initial 
phase-space density but with different realisations (say, by 
choosing a different random-number-generator seed) would 
produce significantly different results when measuring the 
test particles distribution at some final time t. The distri- 
bution of the massive particles, however, would be almost 
statistically identical. This behaviour was observed when we 
first tried to measure Ki(r) using the above method. 

To overcome this difficulty we spread the test particles 
over large 5D regions in phase-space in which /o(r, v) is con- 
stant. We denote these regions by Vo- The assumption here 
is that in such regions the diversity in the test particles tra- 
jectories is large enough to produce stable statistical results. 
This of course has to be checked numerically by running the 
same simulation with a different realisation of the massive 
particles. 

Finally, a similar situation exists in the final snapshot, 
when one wishes to recover the phase-space density of the 
test particles. A single-point measurement tends to fluctu- 
ate both in time and as a function of the initial realisation 
due to Poissonian noise. These fluctuations can be removed 
by measuring the average phase-space density on surfaces 
over which we assume from symmetry considerations that 
the exact phase-space density (i.e., the exact solution of the 
Vlasov equation) is constant. For example, if the equilibrium 
system is spherical and isotropic, we can choose the surface 
to be a sphere around its centre with v = 0. We denote the 
final surfaces by Si. Consequently, the averaged conditional 
probability over the initial region Vo and the final surface 
Si is denoted by (K(Si,V )). 

Even after averaging over Vo and Si our results were 
sensitive to the particular realisation of the massive parti- 
cles that was used. To estimate this sensitivity we measured 
(K(Si, Vo)) from five simulations with different realisations 
of the initial conditions, and calculated the average and scat- 
ter. As discussed in Sec. [3] the typical scatter between the 
different realisations was of the order of 20%-30%. 

We conclude this section by noting that the averaging 
procedure over Vo and Si does not change the prediction 
of the theories. To see why this is the case, let us pass to 
a continuous description of the conditional probabilities by 
replacing the macro-cell coordinate i with the continuous 
coordinate f which specifies the phase-space coordinate of 
the centre of the macro-cell i. Then the conditional proba- 
bility is defined such that K (f , i~)d 6 f is the probability that 
a particle that was initially at r will end up in a small phase- 
space region d 6 f around f . The predictions of the generalised 
NK00 and LB67 theories for the conditional probabilities, 
which are given in Eqs. 1221 1231 , can now be written as 

K {NK) (f,r) = A (JVK) (r)B (iVK) (f) e £ S n , (27) 
K^ LB) (f,r) = A (ifl) (r)B (LB) (f) e c<LB) ^ /oM . (28) 

Then the average conditional probability over Vo and Si is 

(K(Si,V )) = iVor'lSip 1 / d V / d 2 f'K(f',r') 

Jv a J s t 

= | Vol -1 [ d 5 r'K(f',r') . (29) 
Jv 



In the last equality we used the assumption that K[f , t) 
is identical for all f € Si. Plugging the expressions for 
K^ nk ^(t,t) and K^ lb \t,t) into the above equation, and 
using the fact that /o(t) is constant over Vo, we get 

(A, Vb)) 

= (a (jvx) (t)) () B^^i^e^T 11 (30) 

^^ LB \r)) o B (LB \f)e c<LB) ^ f ° , (31) 

with (-) denoting an average over the initial region Vo. 

We thus see that the structure of the conditional prob- 
abilities in both theories remains unchanged. 

3.2 Numerically estimating the phase-space 
density of the test particles 

Estimating the 6D phase-space density of few tens of thou- 
sands particles is by no means a trivial task. A simple box- 
counting procedure with equal volume boxes is impracti- 
cal as the number of boxes in the 6D phase-space over- 
whelmingly exceeds the number of available particles, un- 
less we choose the boxes to be so large that the resul- 
tant resolution is extremely poor. The solution is to use 
an adaptive technique such as the kernel-based technique 
in SPH simulations. In this work, we used the Delaunay 
Tessellation Field Estimator (DTFE) method w hich was in- 
troduced bv lSchaap fc van de Weveaert I (120001) to estimate 
real sp ace densities, and was later adapted bv lArad et al. I 
(2004) for the estimation of phase-space densities. To cal- 
culate the average phase-space density on the Si sphere we 
used a Monte-Carlo integration technique by randomly pick- 
ing Nmc = 5, 000 points on Si and c alculating th e aver - 
age phase-space density at these points. lArad et al. I J2004T) 
demonstrated that f DTFE / f exact ; s approximately given by 
a log-normal distribution, and therefore we calculated the 
average K using a geometrical mean: 

1 JVjwc 

log (if (Si, Vo)) = jj— Yl log[/f TFB (r J )] . (32) 

3=1 

In the above formula 71, ... , tn mc are the Monte-Carlo sam- 
pling points on Si and ff TFE (Tj) is the DTFE estimate for 
the phase-space density of the test particles at Tj - test par- 
ticles that were placed initially in the region Vo. 

To estimate the internal DTFE measurement error 
we used the same Monte-Carlo technique to estimate the 
phase-space density in a synthetic distribution realised using 
100, 000 particles (the same number of test particles as was 
used in the simulations, see next section). The /(r, v) of the 
synthetic distribution was identical to the initial distribution 
of the massive particle in the experiment (see next section) . 
This is a Hernquist-like distribution with Gaussian veloc- 
ity distribution, described by Eqs. (1331361 . We estimated 
the phase-space density on spheres with 0.1 < r < 1.8 and 
v = 0. The results are presented in Fig. We see that 
the DTFE method produces a Poissonian scatter of about 
20%-30% around its mean, which in turn matches the exact 
phase-space density within an error of no more than 10%. 
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Figure 1. The DTFE-measured average phase-space density of 
the massive particles on surfaces where the exact phase-space den- 
sity is constant. Each average phase-space density was calculated 
using 5,000 random points on the surface. The error-bars denote 
the Poissonian error which is typically 30%. However, the error 
in the average density itself is much smaller, typically about 7% 
and it never exceeds 10%. 




Figure 2. An illustration of our initial conditions showing the 
four haloes arranged in a symmetrical tetrahedron. The figure 
depicts the massive particles embedded with the massless test 
particles in the centre of each halo. 



We notice that the average DTFE values are systemat- 
ically smaller than the exact density. This can be explained 
by the fact that we are measuring the phase-space density 
on a surface where v = 0. According to Eq. 13311 in the next 
section, such a surface corresponds to a local maxima in 
/o(r,v) and since the density in a given point is calculated 
by linearly interpolating the phase-space density of nearby 
particles, we expect it to be lower than the exact value. In- 
deed for smaller radii, where the mean separation between 
particles is smaller, the systematic deviation is smaller. Such 
systematic errors, however, will be largely cancelled out from 
the ratio in Eq. . 



4 NUMERICAL SIMULATIONS 

We used a large set of numerical simulations to measure 
(K(Si, Vo)) for 12 initial regions Vb and 5 final surfaces Si. 
All simulations used the same initial phase-space density 
/o(r,v) which was realised using 120,000 massive particles. 
Each of the initial Vo regions was sampled using 100,000 test 
particles. The simulations were run in physical coordinates 
with no cosmological expansion and in dimensionless mode 
with the gravitational constant, unitlength, unitmass and 
unitvelocity all set to unity. 

As was mentioned in Sec. 13.11 we used 5 different re- 
alisation of essentially the same simulation in order to esti- 
mate the sensitivity of the final {K(Si, Vb)) to a particular 
realisation. Therefore, for each Vb we run a set of 5 dif- 
ferent simulations in which we used a different realisation of 
/o(r, v) (by using a different random-number-generator seed 
in the routine that set the initial positions and velocities of 
the massive particles). The initial positions and velocities of 
the test particles in Vb were always the same. Our results, 
which are presented in Sec. |K| are based on the average of 
these different realisations. 

In the following sections we give a detailed description 
of the initial conditions and of the numerical simulations 
themselves. 



4.1 Initial conditions 

4-1.1 Massive particles 

Initially, the system consisted of four identical haloes that 
collapse into each other, producing a single virialised halo. 
Each halo was realised using 30,000 massive particles, giving 
a total of 120,000 massive particles in each simulation. The 
haloes were placed in a symmetrical tetrahedron around the 
origin with their respective centres separated by six length 
units. This setup was chosen to maximise the level of vio- 
lent relaxation and to ensure that there was no preferred 
direction in the final merger of the haloes. The test particles 
were then placed in equi-density regions inside each halo, as 
explained in the next section. A schematic picture of this 
setup is shown in Fig. [5] the four haloes are represented by 
the large spheres while the test particles are represented by 
the smaller spheres which are embedded inside them. 

The haloes were set to be spherical and isotropic with 
Gaussian velocities: 



/o(r,v) = (2^r 3/2 



<j 3 (r) 



We used a Hernquist-like density profile 
p(r) 



r (l + r )3 



(33) 



(34) 



whil e the velocity dispersion a(r ) was set by the Jeans equa- 
tion iBinnev fc Tremaine Il987l) with a vanishing anisotropy 
parameter /3 



dipt??) 
dr 



dr 



This yielded 



2/ \ G 

a ( r = 7 \ ; 

P(r) J r 



M(r')p(r') ^ 



(35) 



(36) 



To increase the amount of violent relaxation in the simula- 
tions, we set the gravitational constant to G = 0.5 in this 
calculation instead of the G = 1.0 that was used in the sim- 
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vo (TO ) 


fo 


VnfO 30) 


0. 


495881 


^0(0.33) 


0. 


416877 


Vb(0.37) 


0. 


337873 


Vo(0.43) 


0. 


258869 


V (0.53) 


0. 


179865 


Vb(0.71) 


0. 


107462 


Vb(0.73) 


0. 


100861 


Vb(0.96) 


0. 


060260 


Vb(1.16) 


0. 


041870 


Vb(1.32) 


0. 


032079 


Vq(1.47) 


0. 


026000 


Vo(1.60) 


0. 


021858 



Table 1. The phase-space density in the 12 initial volumes Vb(ro). 
Each initial volume Vo(ro) is a union of 4 shells around the 4 
haloes. The outer radius of the shell is ro and the inner radius is 
ro/2. The outer radius ro has values in the range 0.30 < ro < 1.60 

ulation. This ensured that each halo would simultaneously 
collapse upon itself while merging with the other haloes. 

4-1.2 Test particles 

We used 12 initial Vo regions. Each region was defined as 
the union of four identical thick shells in real space around 
each one of the 4 haloes. Each shell extended from the radius 
ro/2 to ro, with ro taking 12 possible values from 0.3 to 1.6. 
We denote the Vo region that corresponds to ro by Vo{ro). 
The initial values of ro and the corresponding fo for each Vo 
are given in Table. 14.1.11 

We used 100,000 test particles to sample each Vo region 
by placing 25,000 particles in every shell. The test particles 
were put in locations where the initial phase-space density 
Eq. was exactly fo = (27r)~ 3/2 p(r )/ 'a 3 (r ). The posi- 
tion and velocity of a test particle were chosen in the fol- 
lowing way: first we chose the particle position by randomly 
picking a location within the shell, in such a way that the 
real-space distribution of test particles would be homoge- 
neous and isotropic. Once r was set it implied a value for 
|v| by requiring that the RHS of Eq. 1331 would be equal to 
fo- The directionality of v was then set isotropically. 

Strictly speaking, one should add the contributions to 
/(r, v) from the three other haloes. However, by taking 
ro small enough we have ensured that these contributions 
would never be larger than 10% of fo- 

4.2 The Si surfaces 

The Si surfaces were chosen as two-dimensional shells in real 
space whose centres were at the centre of the equilibrium 
halo with radii of n = 0.3, 0.4, 0.5, 0.6, 0.7. In accordance 
with Sec. 12.41 the measurements were done well within the 
centre of the relaxed halo where the dynamical time is short 
and violent relaxation is believed to be efficient, see Fig. 2] 
and Sec. 14.31 for more details. To specify a particular Si 
surface we use the notation Si(ri). We define the centre of 
the halo by calculating its centre of mass using an iterative 
centre of mass (COM) approach. 

The velocity coordinates of the Si surfaces were chosen 
to be the velocity of the centre of mass, which amounts to 



v = in the centre of mass frame. The underlying assump- 
tion here is that the equilibrium phase-space density of the 
test particles is isotropic when measured in that frame, and 
therefore for f = (r,v) the un-averaged (i.e., local) condi- 
tional probability is given by 

A-(r,r)=A-(|r|,|v|,(r.v),r) . (37) 

Such a function is constant over surfaces with |r| = ri and 
v = 0. Support for this assumption is found in the Poisso- 
nian errors in {K(Si, Vo)) that we measured using the DTFE 
and Eq. 1321 . In all the measurements we obtained an er- 
ror < 30% - which is the typical Poissonian error that was 
found when measuring the phase-space density in the spher- 
ical and isotropic Hernquist-like distribution in Eqs. (1331361 
(sec Fig.Qand Sec. 13.21 . If the phase-space density on these 
surfaces was not constant we would have obtained much 
larger errors. 



4.3 iV-body simulations and numerical effects 

All simulations were run using the par allel version of 
the p ublicly available tree-code GADGET dSpringel et al. I 
2001). The simulations were run on COSMOS, a shared- 
memory Altix 3700 with 152 1.3-GHz Itanium2 processors 
hosted at the Department of Applied Mathematics and The- 
oretical Physics (Cambridge). 

The accuracy of a collisionless GADGET simulation can 
be determined by three parameters: the gravitational soften- 
ing e, the internal time-step accuracy r\ and the cell-opening 
accuracy parameter of the force calculation a. After a num- 
ber of test simulations we settled for the following values 
for the simulation parameters. We used a physically fixed 
force softening of e = 0.02 for both the massive and mass- 
less test particles in the simulations. The half-mass radius 
r m/2 °f the final collapsed halo was ~ 2.7 length units and 
hence the employed softening was e = 7 x 10 -3 r m / 2 - We 
chose timesteps according to At = y / 2r/e/|a| (time-step cri- 
terion of GADGET). Here a is the acceleration and r\ is 
the time-step accuracy parameter which we set to r\ = 0.01. 
In the force calculation we employed the new GADGET cell 
opening criter ion with a high force accuracy of a = 0.001, 
for details see ISprineel et aTlfcOOll) . 

The simulations were ran for t = 70 time units, with 
a typical run taking ~ 70, 000 timesteps to reach the final 
time. In the final snapshot, at least 70% of the test particles 
were found within the half-mass radius, the fraction rising 
to above 90% for runs with low values of ro. At the half- 
mass radius the crossing time was t cross ~ 10, thus ensuring 
that the test particles that were used for the measurement 
were fully relaxed. Further support for this conclusion is 
presented in Fig- 01 which shows the virial ratio —U/2K for 
one of the simulations. Here U is the total potential energy 
and K is the total kinetic energy of the system. For a fully 
relaxed system this ratio should be equal to one (the virial 
theorem). We see that the system fluctuates strongly un- 
til t ~ 20, after which the gravitational potential becomes 
constant and the evolution of the phase-space density (both 
of the massive and test particles) is through phase mixing 
only. Figure0]shows the radial density profile of the massive 
particles at t = 70. The dashed line is the function 4po(r) 
with po(r) being the initial Hernquist-like profile of the four 



A numerical comparison of theories of violent relaxation 9 



1.5 



0.5 - 




20 



40 



GO 



Figure 3. The virial ratio —U/2K as a function of time for one 
of the simulations. U is the total potential energy and K is the 
total kinetic energy. Other simulations give very similar results. 
The dashed line denotes a ratio of 1 - for a fully relaxed system. 
For t > 30 the virial ratio is always between 0.95 and 1.05, which 
are denoted by the dotted lines. The principal fluctuations occur 
at t < 20, whereas at later times the gravitational potential is 
essentially constant and evolution proceeds through mixing. 
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Figure 4. The radial density profile of the relaxed system (solid 
line) at t = 70. The dashed line is the initial radial density profile 
of the 4 haloes scaled by a factor of 4. The good fit between the 
two is in agreement w i th oth er numerical experiments such as 
iBovlan-Kolchin fc Ma I fc004t) . The two arrows denote the rj = 
0.3 and r\ = 0.7 radii which are the minimum and maximum 
radii at which the phase-space density of the test particles was 
measured. 



haloes given by Eq. 1341 . We see that the profile of the re- 
laxed halo is well fitted by the profile of the initial ha loes. 
This is in agreement with Bovlan-Kolchin & Ma ( 2004) who 
found the same behaviour in a head-on collision of two cuspy 
haloes. 

An important test in assessing the reliability of a colli- 
sionless simulation is to calculate to what extent the simu- 
lation is affected by two-body relaxation. There have been 
a number of recent studies on what regions c an be consid- 
ered r el iable against two - body relaxation [i.e. IPower et al. I 
J2003ft: iDiemand et al 1 (120041) 1. We chose here to follow 
iBovl an- Kolchin fc Ma I d2004l) who construct the local two- 



body relaxation time t T as defined in terms of the circular 
velocity Vdrc ~ V c ; rc (r), period T(r) = 2nr/V C irc and the 
number of particles N (or equivalently mass M) interior to 
a radius r as 

t r (r) = T(r vi 



8 111 iV Tyir Vcircij^vir^) 



N 



4 lniVU GM (r 



(38) 



To minimise the effects of two-body relaxation we require 
that t r (r) is longer than the overall simulation time for all 
values of r for which we measure the conditional probabil- 
ities. We find for the range r = 0.3 — 0.7 for which the 
measurement of the final phase-space density is done that 
the relaxation time is t r — 90 — 470, thus concluding that 
the simulation is not affected by two-body relaxation in our 
region of interest. 



5 RESULTS 

Figure shows log (K[Si(ri), Vb(fo)]) as a function of n 
for 6 out of the 12 initial conditions, at two different times: 
ti = 40 time units and ti = 70 time units. To simplify the 
notation, we denote it simply by K(r\,ro). 

K(ri,ro) is the result of averaging both over the 5i(ri) 
surfaces and over the 5 different simulations that were run 
for each initial region V6(ro). We first calculated the Si aver- 
age and then calculated the average of the different realisa- 
tions. Both averages were done using logic". In each realisa- 
tion about 100-700 particles contributed to the Si(ri) aver- 
age. These particles defined the Delaunay cells which were 
intersected by the 5D Si(ri) surface. The errors in Fig. 
are the scatter errors among the different realisations. They 
are typically in the range 20% — 30%, and in general the 
two snapshots agree within that range. We also notice that 
for the lowest values of ro the agreement between the two 
snapshots is better than for higher values of ro. This is an 
indication that the phase-space structure of the low ro test 
particles is more relaxed than that of the higher ro. In addi- 
tion, for lower values of ro, the spatial distribution is more 
concentrated, having more test particles with smaller orbits. 
Such particles have shorter dynamical times and thus tend 
to relax more quickly than particles on outer orbits. It is 
therefore not surprising that K(ri,ro) with lower ro are less 
fluctuating than those with higher values of ro- 

Figure [(J shows the log of the ratio R(ri,rj) = 
K(ri,ro)/ K(rj ,ro) calculated from the average K(r, ro) at 
t = 70. This is the main result of this paper. The upper 
plots show log R as a function of fo for comparison with the 
LB67 theory, whereas the lower plots show the same log R 
as a function of l//o for comparison with the NK00 theory. 
In both groups the pairs (ri,rj) are (0.3,0.5), (0.5,0.7) and 
(0.3, 0.7). Other combinations of radii give similar results. 

In the upper plot we see that log i?(r;, r 3 ) increases until 
/o ~ 0.25 after which it becomes approximately constant. 



On the other hand, in the lower plot log R(rt 



decreases 



until l//o ~ 25, where again it saturates. Putting it all 
together we see that log i?(r;, r 3 ■) is approximately constant 
for fo < 0.04, then it increases (as a function of fo) until 
fo ~ 0.25 and then it saturates once more, log R(n,rf) is 
therefore highly non-linear both as a function of fo and as 
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a function of l//o- As was discussed in Sec. 12.41 a linear 
behaviour is a necessary condition for either theories to be 
correct. 

To quantify the above departure from non-linearity, 
we fitted a straight line for each of the size plots using a 
least-mean-square procedure. As expected, the results, are 
decisively against both theories. The quality of the fits is 
very bad, as can be clearly seen from the plots. The re- 
duced x 2 °f the fits for the Lynden-Bell theory are (left 
to right, top of Figure 01 Xlb = 1-31,2.51,6.59 and for 
the Nakamura theory (left to right, bottom of Figure |SJ 
Xnk = 4-06, 2.43, 10.00. The lower X 2 for the first two fits is 
due to the small differences between and rj , which causes 
Rij to be very close to unity and consequently log Rij to ap- 
proach zero, and thus be better approximated by a straight 
line. On the other hand, in the third pair where the differ- 
ence between r» and rj is maximal, Rij spans a wider range 
of values, and consequently both theories produce a very 
poor fit of x\b = 6.59 and Xnk = 10.00. This clearly shows 
that both theories fail miserably the test of linearity of log R 
in f (LB67) and l// (NK00). 



6 SUMMARY AND CONCLUSIONS 

We have performed a series of iV-body simulations to test 
the validity of the LB67 and NK00 theories. Unfortunately, 
due to a limited amount of time and computer resources we 
did not perform simulations with different /o(r, v) or with 
higher number of particles. Different initial conditions and 
different resolutions may yiel d qualitatively different equili b- 
rium states [see, for example. iMerrall fc Henriksen I J2003)]. 
The running of additional simulations is left for possible fu- 
ture work. However, in light of the strong non-linearity of 
the plots in Fig. |S| and the use of a very weak and general 
condition to test the theories, we believe that our conclu- 
sions will not be altered by additional simulations. 

The main aspects in which our numerical experiment 
differ from previous studies are: 

• We used 3D simulations whereas previous studies con- 
centrated on ID simulations. 

• The experiment was explicitly designed to check which 
one of the two possible formulas of entropy is more correct: 
the LB67 formula which is derived using an equal-volume 
discretisation or the NK00 formula which is derived using 
an equal-mass discretisation. 

• To distinguish between the two theories we used a very 
weak condition on the conditional probabilities Ki(t), which 
was estimated by measuring the phase-space densities of sets 
of test particles. Therefore we did not need the full analyt- 
ical solution of the theories with all its computational and 
conceptual difficulties. 

• Previous attempts to verify the LB67 theory were done 
with the water-bucket initial conditions, in which the initial 
phase-space density has only one level. Here, as we needed to 
distinguish between the LB67 and NK00 theories, the initial 
phase-space density covered a continuous range. 

The results of the experiment are summarised in Fig. |S| 
They provide very strong evidence against the LB67 and 
NK00 theories. As was discussed in Sec. 12.41 the linearity of 
logi? in /o or l//o is a very basic and weak requirement of 



both theories. Therefore the non-linearity of all the plots in 
Fig. |S| must be attributed to the failure of the most basic 
assumption in these theories, i.e., the maximisation of en- 
tropy. Indeed, we cannot explain the failure of the theories 
by the existence of additional conserved quantities that de- 
pend on the /(r, v) (such as the total angular momentum 
vector L) since we have used the final /(r, v) itself as a con- 
straint. We also cannot argue that our measurements reflect 
the incomplete relaxation of the outer parts of the system 
since they are done in the innermost parts of the system 
(less than l/2r m / 2 ) and in addition the linearity condition 
is independent of the full solution of the theories that as- 
sumes complete relaxation in all parts of the system. As was 
explained in Sec. 12.41 the weaker the condition is, the deeper 
is the failure of the theories if they do not pass it - and this 
is the bottom line. 

In some respect this is a disappointing result as it shows 
that neither of the theories is even remotely correct. Ad- 
ditionally we are unable to decide which definition of the 
entropy is the "right" one as they both seem to perform 
equally bad. 

On the other hand, one may find some sort of comfort 
in the fact that both theories fail, as there is no solid a priori 
theoretical argument against either equal- volume discretisa- 
tion or equal-mass discretisation, and it is not obvious why 
they should produce such different re sults in the first place- 
Addit ionally, as was shown recently bv lArad fc Lvnden-Bell I 
(2005), both theories contain some sort of self-inconsistency 
since they are both non-transitive. Knowing that the theo- 
ries are fundamentally wrong empirically thus solves many 
of these problems or at least makes them less relevant. 
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Figure 5. The average log K(r, ro) for 6 out of the 12 possible values of the initial radius ro. The dashed line is the average at t = 40 
whereas the smooth line is the average at t = 70. The average was calculated using 5 different runs of the same simulation with different 
realisations of the initial conditions. The statistical error is typically smaller than 20%. The agreement between the t = 40 and t = 70 
snapshots is also typically within an error of 20%. It is particularly good for the lowest values of ro, indicating that these test particles 
are more relaxed that the ones with larger values of ro- 




-0.2 




Figure 6. The logarithm of the ratio R(ri, rj) = K(rt, ro)/K(rj , ro) as a function of /o and of l//o, with /o being the initial phase-space 
density of the test particles. The K(rj, ro) that was used to calculate the ratio is the average of 5 simulations, and is shown in Fig. 151 
For all plots a straight line was fitted using a least-mean-square procedure. The resultant x 2 for the three upper plots (LB67) is (left to 
right) x 2 = 1-31,2.51,6.59 and for the lower plots (NK00) \ 2 = 4.06,2.43,10.00. All plots show a very strong non-linear behaviour, in 
contradiction with both LB67 and NK00 theories. 
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